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Abstract 

(N 

^ ' While kinetic Monte Carlo simulations can provide long-time simulations of the 

dynamics of physical and chemical systems, it is not yet possible in general to iden- 
■^ i tify the inverse Monte Carlo attempt frequency with a physical timescale in any 

^^ I but the simplest systems. Here we demonstrate such an identification by comparing 

—1. ■ simulations with experimental data. Using a dynamic lattice-gas model for the elec- 

^^ , trosorption of Br on Ag(lOO), we measure the scan-rate dependence of the separation 

between positive- and negative-going peaks in cyclic voltammetry and compare sim- 
ulated and experimental peak separations. By adjusting the Monte Carlo attempt 
frequency, good agreement between simulated and experimental peak separations is 
'^ I achieved. It is also found that the uniqueness of the determination depends on the 

^ ' relative values of the adsorption/desorption and diffusion free-energy barriers. 

O 
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1 Introduction 



At present, kinetic Monte Carlo (KMC) simulation is virtually the only com- 
putational method that enables numerical study of the dynamics of physical 
and chemical systems on macroscopically relevant timescales, anywhere from 
microseconds to millions of years [1,2,3,4,5]. However, the method is essentially 
a stochastic approximation to underlying classical or quantum mechanical pro- 
cesses on finer space and time scales [6,7]. It thus suffers from the problem 
that the basic MC timescale is often difficult to relate to an underlying phys- 
ical timescale. While such a timescale could, in principle, be calculated from 
comparison with calculations at these finer scales [6], this is in practice rarely 
possible for electrochemical systems. Due to the complexity of the interactions 
with the solution, ab-initio methods can construct a reasonable horizontal 
corrugation potential [7] but give limited knowledge about the shape of the 
potential in the direction perpendicular to the surface. Moreover, water has ef- 
fects both in terms of damping and in the shape of the adsorption/desorption 
free-energy barrier, which probably corresponds to a reconstruction of the sol- 
vation shell. While analytic macroscopic theories can be derived in terms of 
MC parameters [8], such theories cannot directly predict the values of those 
parameters. A rough estimate of the overall reaction rate constant can be ob- 
tained by applying standard electrochemical techniques [9] to the experimental 
data. Yet, the overall reaction rate gives little information, if any, about the 
free-energy barrier heights for the different processes of adsorption/desorption 
and surface diffusion. 

Here we present an alternative approach: comparison of KMC results with 
time-dependent experimental results. In particular, we compare KMC simula- 
tions of a lattice-gas model with experimental results for cyclic-voltammetry 
(CV) studies of the electrosorption of Br on single- crystal Ag(lOO) surfaces. 
The lattice-gas model represents the long-lived configurations of adsorbed Br, 
and an MC step corresponds to an attempt at hopping across a saddle point 
in the free-energy landscape to a new configuration [7]. 

The Br/Ag(100) system has a phase transition (in the two-dimensional Ising 
universality class) between a disordered phase at more negative potentials and 
an ordered c(2 x 2) phase at more positive potentials [10]. The phase transition 
is associated with a divergence of the coverage fiuctuations, corresponding to 
a peak in the cyclic voltammogram. The same phase transition has also been 
observed for Cl/Ag(100) in electrochemical [11] and ultra-high cavuum (UHV) 
environments [12,13]. 

Recent static and kinetic MC studies have been used to investigate the phase 
ordering and disordering mechanisms in cyclic-voltammetry (CV) and sudden 
potential-step experiments for halide adsorption on Ag(lOO) [14,15]. As the 



CV scan rate is increased, the system is driven further from equihbrium. As 
a result, there is a widening of the separation between the positions of the 
CV peak for the positive-going and the negative-going scan of the electrode 
potential (in experiments) or the electrochemical potential (in simulations). 
Experimental peak separations have been measured for different sweep rates 
of the electrode potential [16]. When these experimental values, with time 
measured in seconds, are compared to the simulations, with time measured in 
Monte Carlo steps per site (MCSS), a physical time can be associated with 
the inverse MC attempt frequency. This was previously attempted by Mitchell 
et al. [17], but simulations at that time could not achieve peak separations 
within the experimental range. Due to the increase in computer power and to a 
new mean-field enhanced simulation method that we developed for calculating 
long-range interactions [11], it is now possible to simulate peak separations well 
within the experimental range. 



2 Lattice-gas Model 



We employ a LG model similar to that used by Koper [18,19] and Mitchell, 
et al. [14,15]. The Br ions adsorb on four-fold hollow sites of the Ag(lOO) sur- 
face [7,20]. The model is defined by the grand-canonical effective Hamiltonian, 

N 

Ti = - ^ (pijCiCj -Jl^Ci , (1) 

where I]i<j is a sum over all pairs of sites, (pij are the lateral interactions 
between particles on the ith. and jth lattice sites measured in meV/pair, /I is 
the electrochemical potential measured in meV/particle, and A^ = L^ is the 
total number of lattice sites. The local occupation variables q can take the 
values 1 or 0, depending on whether site i is occupied by an ion (1) or empty 
(0). 

The simulations were performed on LxL square lattices, using periodic bound- 
ary conditions to reduce finite-size effects. The interaction constants (pij be- 
tween ions on sites i and j a distance Tj^ apart (measured in units of the 
Ag(lOO) lattice spacing, a = 2.889 A [10]) are given by 

(pij = -CX) X dr.^^i H 3 , (2) 

where 5^. ,i is a Kronecker delta with the infinitely negative value for Vij = 1 
indicating nearest-neighbor exclusion. Negative values of 0nnn denote long- 
range repulsion. The interactions for large Vij are most likely of electrostatic 



dipole-dipole nature [21], but may also have a component mediated by the Ag 
substrate [22] that may not be uniformly repulsive. 

Assuming a sufficient concentration of counterions in the electrolyte, the elec- 
trochemical potential JI is, in the dilute-solution approximation, related to the 
bulk ionic concentration C and the electrode potential E (measured in mV) 
as 

C 

JI = JIq + k^T In — - e'-fE , (3) 



where /Iq is an arbitrary constant, Cq is a reference concentration (here 1 mM), 
k-oT is Boltzmann's constant times the temperature, and e is the elemen- 
tary charge unit. The electrosorption valency 7 [23,24] (here assumed con- 
stant [11,21]) corresponds to the fraction of the ionic charge transferred dur- 
ing the adsorption process. We use the sign convention that 7^ > favors 
adsorption. 

N 

The adsorbate coverage is defined as ^ = A^~^ X/*^*- ^^ '^^'^ ^^ experimentally 

obtained by standard electrochemical methods, as well as from the integer- 
order peaks in surface X-ray scattering (SXS) data [10,25]. The derivative of 
the coverage with respect to the electrochemical potential, d^/d/I, is propor- 
tional to the current density in an experimental CV as described in section 4. 

We have previously estimated the next-nearest-neighbor lateral interaction en- 
ergy (0nnn ~ "21 meV) and the electrosorption valency (7 ~ —0.72), based on 
fits of simulated equilibrium adsorption isotherms to experimental data [11]. 
Remarkably, despite the very different environments this value of 0nnn, esti- 
mated from the shape of the adsorption isotherms, is consistent to within 10% 
with the one obtained from the temperature dependence of the critical cover- 
age for Cl/Ag(100) in UHV [12,13]. In the next section we detail the kinetic 
MC approach. 



3 Kinetic Monte Carlo Method 



Kinetic Monte Carlo simulations were performed on systems with L = 128, 
and the absence of significant finite-size effects was verified by additional sim- 
ulations for L = 64 and 256. The kinetic MC simulation proceeds as follows. 
We randomly select a lattice site i, and depending on the occupation q, dif- 
ferent moves are attempted. If q = 0, only adsorption is attempted, while if 
Cj = 1, 9 different moves are proposed: desorption, diffusion to each of the 
4 nearest-neighbor sites, and diffusion to each of the 4 next-nearest-neighbor 



sites. Next, a weighted list for accepting each of these moves is constructed us- 
ing Eq. (5) below, to calculate the probabilities i?(F|I) of the individual moves 
between the initial state I and final state F. The probability for the system to 
stay in the initial configuration is consequently -R(I|I) = l—TipjLiR(F\l) [15]. 

Using a thermally activated, stochastic barrier-hopping picture, the energy of 
the transition state for a microscopic change from an initial state I to a final 
state F is approximated by the symmetric Butler- Volmer formula [2,26,27] 

U,, = ^l±^ + A,, (4) 



where Uj and Up are the energies of the initial and final states, respectively, 
Ta is the transition state for process A, and Aa is a "bare" barrier associated 
with process A, which here can be one of nearest-neighbor diffusion (Ann), 
next-nearest-neighbor diffusion (Annn), or adsorption/desorption (Aa/d)- 

The probability for a particle to make a transition from state I to state F is 
approximated by the one-step Arrhenius rate [2,26,27] 

^,p|,)..e.p(-^)exp(-i|^). (5) 



where u is the attempt frequency, which sets the overall timescale for the simu- 
lation. This formalism assumes that all processes A have the same MC attempt 
frequency of 1 MCSS~^. Although this assumption may not be true in gen- 
eral [28], because different processes involve different potential wells and thus 
different vibrational frequencies, it can be used to find an effective timescale 
conversion between MC time and physical time, given the MC attempt fre- 
quencies. The MC time unit, one MCSS, corresponds to the attempted update 
of A^ randomly chosen sites. For a simulation with scan rate p (in meV/MCSS), 
JI is incremented by (p x 1 MCSS) meV every MCSS until it reaches its final 
value, and then decremented back to its initial value. 

Since the value of Up — Ui determines 7^ in Eq. (5), the approximations made to 
calculate the large-r contributions to the pair sum in Eq. (1) are important. 
In our simulations, to calculate the energy changes we included the exact 
contributions for particle separations up to r^j = 3, while using a mean-field 
approximation for larger separations [11]. ^ 



^ Omission of the mean-field part, which could be viewed as a simple-minded way 
to introduce electrostatic screening, would increase the value of (/'nnn that best fits 
the equilibrium isotherms by 5 — 10% [11,21]. 



4 Results 



Starting at an initial potential of /I = —200 meV and room temperature 
{P = l/k-oT = 0.04 meV~^), the MC process discussed in Sec. 3 was repeated 
until JI reached +600 meV, and then decremented at the same rate back to 
—200 meV. The coverage isotherms were computed using scan rates ranging 
over four decades, from p = 3 x 10~^ to 0.1 meV/MCSS, and for Aa/d = 150, 
175, 200, 250, 300, 350, and 400 meV. The other barriers, An„n = 200 meV and 
Ann = 100 meV, were kept constant at the values determined by comparison 
with density-functional theory calculations [7,20]. The hysteresis loops for 6 
as a function of E are shown for Aa/d = 300 meV in Fig. 1 for different scan 
rates. Each loop was averaged over eight independent simulation runs. 

Next, the Savitzky-Golay method [29,30] with a second-order polynomial and 
a window of 51 points was used to obtain the smoothed numerical derivative 
dO/dp, which is proportional to the simulated CV current density j [15]: 

, _ 7^62 dedE 



where Ag is the area of one adsorption site. Consequently, the differential 
adsorption capacitance per unit area becomes 

C = ^ = 2!^^ (7) 

dE/dt A, dji' ^ ^ 



As seen in Fig. 2, the difference between the positions of the positive-going 
and negative-going peaks depends on the scan rate. This difference is due to 
two physical effects. First, jamming due to the nearest- neighbor exclusion, 
which is only slowly alleviated by lateral diffusion [15], causes the coverage to 
lag behind its equilibrium value. Second, critical slowing down also causes the 
coverage to equilibrate slowly near the phase transition at the CV peak [15]. ^ 

By fitting the separation of the positive-going scan peak {Ep) and the negative- 
going scan peak (En) of the simulated CVs to experimental data for Ep — En 
vs dE/dt [16], we can extract a physical time r corresponding to 1 MCSS 
(the inverse MC attempt frequency, r = 1/z/) for each value of the adsorp- 
tion/desorption free-energy barrier height Aa/d- A x^ which measures the 
square of the horizontal differences between the experimental curves and the 



^ If j is plotted vs 6, rather than vs E or JI, all the peaks fall in the range 0.37 < 
6 < 0.40 for positive-going scans and 0.33 < 9 < 0.37 for negative-going scans. 
This is consistent with the values of the hard-square jamming coverage and critical 
coverage [12,13,15]. 



simulated peak separations, multiplied by trial values of r, was calculated. The 
value of T which minimized x^ 'was taken to be the best-fit value. See Fig. 3. 
For Aa/d < 350 meV, the simulations were fit to the experimental data points, 
while for Aa/d = 400 meV the simulations were fit to the best-fit simulated 
data points for Aa/d = 300 meV since the peak separations attainable for 
Aa/d = 400 meV are not within the experimental data range. The fits for 
most values of Aa/d coincide to within the accuracy of the experimental data 
and our statistics. Yet, focusing on the experimental data range (inset in 
Fig. 3), suggests that Aa/d = 175 meV and r = 5.3 x 10~^ s fits best to 
the experimental data. This distinction relies mainly on the two experimental 
data points corresponding to the highest scan rates and thus would be better 
founded if there were more experimental data in that range. In addition, since 
the simulated curves for Aa/d > 200 meV practically coincide, it is only for 
Aa/d < 200 meV that a clear distinction can be made among the different 
Aa/d values, even with more accurate experimental data. 

For Aa/d > 200 meV, the physical time r corresponding to one MCSS is 
simply related to Aa/d as seen in Fig. 4. The relationship between r and Aa/d 
is related to the Arrhenius form: 

r = roexp(-/5Aa/d) , (8) 



or 



logio (r/s) = ln(ro/s) - /3Aa/d / In 10. (9) 



Plotting log^Q (r/s) vs Aa/d results in a straight line with a slope of 
—0.0388 meV~^/lnlO (excluding Aa/d = 150 and 175 meV), in very good 
agreement with what is expected from Eq. (9) with the inverse temperature 
used, (3 = 0.04 meV""*^. We also find that, as Aa/d decreases and becomes 
comparable to Annn and Ann, diffusion becomes relatively more important in 
determining the overall timescale, and the dependence of r on Aa/d deviates 
from Eq. (9). This deviation makes differentiation among different values of 
Aa/d easier because in this limit a change in Aa/d cannot be compensated by 
a change in r, and the details of the dynamics directly influence the peak 
separation. 



5 Conclusions 



By comparing with experiments simulations at potential-scan rates sufficiently 
slow to produce peak separations that fall within the experimental range, we 



were able to extract a physical timescale associated with the inverse MC at- 
tempt frequency r = 5.3 x 10~® s, a value much larger than normally expected. 
This may be the result of relatively flat potential minima for the adsorption 
process, possibly related to reorganization of the ion hydration shells. An- 
other possible reason is the assumption that all processes in the MC have the 
same attempt frequency; different processes could, in general, have different 
MC attempt frequencies. Thus r would correspond to an effective inverse MC 
attempt frequency for an overall process. 

For values of the adsorption/desorption free-energy barrier Aa/d that are large 
compared to the diffusion barriers (Aa/d > 200 meV), the relationship between 
r and Aa/d is consistent with the Arrhenius law. Under such conditions the 
process of Br adsorption/desorption controls the dynamics in this electrochem- 
ical system, and the possible difference in attempt frequencies for different 
processes becomes less important. For situations in which Aa/d is comparable 
to the values of the lateral diffusion of adsorbates on the substrate, deviations 
from the Arrhenius law appear. Then no process is dominant. In conclusion 
we have here shown that simulations which measure the competition between 
these processes can be used to distinguish between different Aa/d values by 
comparison to experiments. 
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Fig. 1. The coverage as a function of the electrode potential E for Awd = 300 nieV 
and different scan rates. 
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Fig. 2. Simulated adsorption capacitance cyclic voltanimograms for Awd = 300 meV 
and different scan rates. 
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Fig. 3. Fits of kinetic MC simulations to experimental peak separations using dif- 
ferent values for the adsorption-desorption barrier, A^/d- Zooming into the exper- 
imental data range (inset) suggests that the best-fit value for A^/d is 175 meV, 
corresponding to r ~ 5.3 x 10~^ s. The lines are guides to the eye. 
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From fits to experiment 
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Fig. 4. Arrhenius behavior of r ~ Toe~^ ^Z'* as a function of Aw^. The slope of 
the hne representing the hnear regression of log^g i'^/^) ^'^^ A^/d > 200 meV is 
(—0.0388/ In 10) nieV"^, in very good agreement to what is expected from the 
Arrhenius equation: /? = 0.04 meV~^. As A^/d decreases closer to the values of 
A„nn and Ann, T deviates from the Arrhenius behavior. 
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